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ABSTRACT 


The  objective  of  this  thesis  was  to  study  whether  simplified  three 

dimensional  finite  element  models  could  be  used  as  a  micromechanical  model  for 

analysis  of  stresses  at  the  fiber  and  matrix  level  while  providing  smeared  composite 

properties  for  the  global  structural  analysis.  The  mathematical  model  was 

programmed  using  the  MATLAB  engineering  software.  For  considering  interface 

cracks  between  the  fiber  and  matrix,  springs  were  modeled  in  the  connecting  nodes 

between  the  fiber  and  matrix  on  each  of  the  three  degrees  of  freedom.  When 

results  were  compared  with  experimental  data,  findings  proved  that  the  simplified 

three  dimensional  finite  element  models  could  simulate  the  material  property  of  the 

graphite/epoxy  composite  very  efficiently.  _ _ _ 
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I.  INTRODUCTION 


Fibrous  composite  material  is  useful  for  structural  application  because  of  its 
beneficial  material  properties  such  as  high  ratios  of  stiffiiess-to-weight  and 
strength-to-weight.  In  fibrous  composite,  fibers  are  embedded  in  a  matrix 
material.  The  matrix,  besides  holding  the  fibers  together,  has  an  important 
function  of  transferring  the  applied  load  among  the  fibers. 

In  order  to  analyze  the  fibrous  composite  structure,  it  is  necessary  to  obtain 
material  properties  for  the  composite.  The  material  properties  can  be  obtained 
from  experiments.  However,  some  material  properties  are  not  easy  to  measure 
directly  from  the  physical  test.  In  addition,  if  there  is  a  different  fiber  volume 
fraction,  the  composite  material  properties  also  vary  even  if  the  same  fiber  and 
matrix  materials  constitute  the  composite.  As  a  result,  it  is  useful  to  predict  the 
material  properties  of  a  composite,  given  the  component  properties  and  their 
geometric  arrangement. 

The  objective  of  this  study  is  to  develop  a  micromechanical  model  to 
determine  the  composite  properties  when  the  fiber  and  matrix  materials  are  known 
as  well  as  the  fiber  volume  fraction.  There  have  been  many  studies  for  this 
purpose  so  far.  However,  most  of  the  studies  were  two-dimensional  models  so 
that  they  provided  limited  information.  A  small  number  of  three-dimensional 
models  were  proposed  [Ref  1-3].  These  three-dimensional  micromechanics 
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models  assumed  perfect  bonding  at  the  fiber/matrix  interface.  In  order  to  model  a 
more  general  case  including  partial  interface  crack,  a  simplified  three-dimensional 
finite  element  model  is  developed  in  this  study  and  the  results  are  compared  to 
available  experimental  results. 

Figure  1  shows  the  micromechanical  modeling  process.  The 
micromechanics  model  considers  a  fiber  surrounded  by  a  matrix  material.  In  the 
following  chapters,  both  analytical  and  finite  element  micromechanical  models  are 
presented,  as  well  as  their  results  for  perfect  interfacial  bonding  and  partial 
interfacial  cracks.  Conclusions  and  recommendation  follow. 
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Figure  1  (a)  to  (d)  shows  macromechanical  analysis  of  a  laminate  composites  and 
micromechanic  analysis  of  a  thin  unidirection  lamina. 
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11.  MICROMECHANICAL  MODELS 


A.  ANALYTICAL  MICROMECHANICAL  MODEL 


Kwon's  micro-mechanics  model  considers  a  unit  cell  consisting  of  the  fiber 
and  surrounding  matrix  material.  The  fiber  is  assumed  to  have  a  square 
cross-section.  Based  on  symmetry  consideration,  a  quarter  of  the  unit  cell  is 
divided  into  four  subcells.  The  size  of  subcell  depends  on  the  fiber  volume 
fraction,  as  shown  in  Figure  2.  The  stress  continuity  at  the  subcell  interfaces  is 
expressed  as 


CT22-CI22?  C>22“^22^  C733-^33i  0^33-033 


<^12 


-  <^12  : 


c  d 

On=^n- 


ai3  =0' 


13  : 


<^13  -  *^13 


(2.1) 


abed 

-  ^23  =  C723  =  C723 

and  the  strain  compatibility  is  assumed  as 


Sii  -  Zii  -  -  Sii 

S22 +£22  =^22 +S22,  833  +  £33  =  £33  +  £33  (2.2) 


S12  +£n  = 


'12 


12  +  £12 : 


'13 


+  £i3  -8*3  +8^3 


The  constitutive  relation  for  each  subcell  is  expressed  by  the  generalized 
Hooke's  law  of  Equation  (2.3) 
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(2.3) 


ij,  k,l=  1 , 2, 3  and  a  =  a,  b,c,d 

The  size  of  each  subcell  depends  on  the  fiber  volume  fraction.  Kwon 
expressed  the  composite  stresses  and  strains  as  a  function  of  the  subcell  stresses, 
the  subcell  strains,  and  the  fiber  volume  fraction  as  shown  in  Equation  (2.4) 

/J=  1,2,3 

where  Oy  and  Sy  are  composite  stresses  and  strains,  o“  and  are  stresses  and 

strains  of  a  subcell  (a  =  a,  b,  c  or  d),  Vf  is  the  fiber  volume  fraction,  and  subcell 

'a'  represents  a  fiber  and  the  rest  of  the  subcells  are  matrix.  The  composite 
stresses  or  strains  are  determined  by  the  volume  average  of  the  stresses  and  strains 
of  the  subcells. 

Instead  of  the  second  and  third  expressions  in  Equation  (2.2),  we  may 
assume  the  following  deformation  compatibility 

JVfsi,  +  (l  -  +  {2.5) 

+  (l  -  e),  =  +  (l  -  4 


7 


This  is  the  primary  difference  from  Kwon's  analytical  model. 

Some  algebraic  operations  are  conducted  for  Equation  (2.1)  through  (2.5)  to 
relate  the  subcell  strains  to  cell  strains  (composite  strains)  as  given  in  Equation 
(2.6). 

M{8“}  =  {s^}  1,2,3  (2.6) 

Inverting  the  coefficient  matrix  [A],  Equation  (2.6)  can  be  also  written  as 

{8“}=(.4r‘w)  (2.7) 

Use  of  Equation  (2.3),  the  first  Equation  of  (2.4)  and  Equation  (2.7)  results  in 
relation  between  composite  stresses  and  strains,  That  is 

{Oy}  =  [C]  {8;y}  (2.8) 

where  matrix  [C]  is  the  smeared  material  property  matrix  for  the  composite. 

Assuming  orthotropic  materials  for  the  fiber  and  matrix,  normal  stress  and 

strain  components  are  uncoupled  from  shear  stress  and  strain  components.  In  this 
case,  the  equations  relating  subcell  strains  to  cell  strains  for  the  normal  strains  are 
given  in  Appendix  A. 

The  comparison  between  analytical  and  experimental  results  is  shown  in 
Figure  3,  for  the  transverse  elastic  modulus  of  a  graphite/epoxy  composite  whose 
material  properties  are  shown  in  Table  1.  The  modified  analytical  solution  agrees 
slightly  better  with  the  experimental  data  than  the  original  analytical  solution  in 
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TRANSVERSE  ELASTIC  MODULUS 


Ref.  [1].  However,  other  properties  such  as  longitudinal  elastic  modulus,  shear 
modulus,  and  Poisson's  ratio  were  not  different  between  the  two  analytical 
solutions  which  agreed  well  with  the  experimental  data. 


Figure  3  Transverse  elastic  modulus  of  a  graphhe/epoxy  composite 
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Table  1  Characteristic  Of  Matrix  And  Fibers  In  The  Graphite/Epoxy 
Composites 


Properties 

Graphite  Fiber 

Epoxy  Matrix 

Longitudinal  elastic 
Modulus  (Gpa) 

232 

5.35 

Transverse  ealstic 

Modulus  (Gpa) 

15 

5.35 

Longitudinal  Poisson's 
Ratio 

0.275 

0.354 

Transverse  Poisson's 

Ratio 

0.49 

0.354 

Inplane  Shear  Modulus 
(Gpa) 

24 

NA 

*  Graphite  fiber  is  transversely  isotropic  material 


Epoxy  matrix  is  isotropic  material 
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B.  FEM  MICROMECHANICAL  MODEL 


The  finite  element  method  [Ref.  4,  5]  is  a  numerical  technique  which  gives 
an  approximate  solution  to  differential  equations  that  model  physical  and 
engineering  problems.  The  finite  element  method  requires  a  problem,  which  is 
defined  in  a  geometrical  space  (or  domain),  to  be  subdivided  into  a  finite  number 
of  smaller  regions  (mesh).  Over  each  finite  element,  the  unknown  variables  are 
approximated  by  using  nodal  variables;  these  functions  can  be  linear  or  higher 
order  polynomial  expressions  that  depend  on  the  node  numbers  and  their  locations 
used  to  define  the  finite  element  shape.  The  governing  differential  equations  are 
transformed  into  a  matrix  expression  over  each  finite  element  and  the  matrix 
expressions  are  summed  ("assembled")  over  the  entire  problem  domain.  As  a 
result,  a  set  of  linear  equations  is  obtained  in  terms  of  the  unknown  nodal 
variables.  In  this  section,  the  three  dimensional  finite  element  method  is  introduced 
followed  by  a  discussion  of  the  programming  procedure  : 

1.  Three  Dimensional  Finite  Element  Method 

The  Galerkin  form  of  the  weighted  residual  procedure  is  adopted  to 
formulate  the  finite  element  method.  The  idea  of  the  weighted  residual  method  is 
that  we  multiply  the  residual  by  a  weighting  function  and  enforce  the  integral  of 
the  weighted  expression  to  vanish.  The  finite  element  formulation  for 
three-dimensional  elasticity  can  be  shown  below  : 
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The  equilibrium  equations  are 


dCx  I  ^xy  8Xxz  A 

dy  ^  dz 

d^yx  CGy  dXyz  „ 

dxrx  ^:y  ,  502  _  A 
dx  ^  dy  ^  dz 

where  the  shear  stresses  satisfy  the  following  relation 


(2.9) 


Xxy  —  Xyx 

Txr  =  Xzc  (2-10) 

'^yz  ~  '^zy 

due  to  symmetric  stress  tensor. 

The  first  step  in  the  FEM  is  to  convert  the  partial  differential  equations  into 
an  integral  expression  by  applying  the  method  of  weighted  residual  to  the  three 
dimensional  stress  equilibrium  equations  (2.9).  Each  equation  is  multiplied  by  a 
weight  function  (i=  1,2,3),  which  is  continuous  over  the  physical  domain  of  the 
problem.  The  resulting  weighted  residual  equations  are 


12 


Applying  integration  by  parts  and  combining  the  three  expressions  yields 


where 


r 


dWi  dWi  dWi 

Ox  dx  +  ^  +  ^XZ  0^ 

dWi  ,  dWz  ,  dlV2 

dx  dy  +  dz 

dw.  dfVi  ,  ^  dWi 

dx  V  dy  dz 


>dv=l< 

ds 


(2.12) 


(j)x  —  Oxt^x  "I"  "^xy^y  txz^z 

^y  =  Txynx  +  CyWy  +  (2. 13) 

^z  =  Xxz«x  +  T>-z«y  +  CTzWz 

are  boundary  tractions. 


Let  'u',  V,  and  V'  represent  displacements  in  the  x,  y  and  z  direction. 


respectively.  Then,  the  relationship  of  strain  to  displacement,  assuming  small 
displacements,  can  be  written  as 


o  —  ^  o  —  —  a  —  — 

~  dx  ~  ~  dz 


(2.14) 


yxy  = 


_  (  dw  ,  dv  ,  .. 

dv  dx  J  ’  yy^  I  5v  5z  j  ’  y^ 


du  .  dw 
dz  dx 


The  constitutive  relation  between  stress  and  strain  is 

{o}=[/)]{8}  (2.15) 

the  material  property  matrix  [D]  for  a  composite  material  is  shown  in  detail  in 
Appendix  B. 
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Equation  (2.12)  can  be  written  in  matrix  form  by  substituting  the 
strain-displacement  Equation  (2.14)  into  the  constitutive  relation  given  by  (2.15). 
The  resultant  expression  is  shown  in  Equation  (2.16). 


li 


fL  0  0 

OX 

0  5  0 

dy 

n  A  ^^3 


dw^ 

dy 

dWi 

dx 

0 


0 

dW2 

dz 

dW-i 

dy 


dWx 

dz 

0 

dW^, 

dz 


dx 

0 

0  ’ 

- 

0 

d 

dy 

0 

/- 

0 

0 

£ 

dz 

0 

[D] 

£ 

< 

dy 

dx 

0 

£ 

V 

dz 

dy 

0 

£ 

dz 

dx 

U 

V 

w 


dV^W, 


(t>xWx 

^yny 

(i)z«z 


ds 


(2.16) 

Eight  noded  isoparameric  rectangular  parallelepiped  elements  are  used  for 
the  formulation.  Each  node  has  three  orthogonal  displacements,  'u',  'v\  and  ’w’,  so 
that  each  element  has  a  total  of  24  degrees  of  freedom  (dof).  In  the  next  step  of  the 
FEM,  we  can  define  the  displacements  in  terms  of  polynomial  shape  fimctions. 


named  //,  ,  as  follows 


u=tHiUi  ■  v=i.HiVi  ;  w=i.HiWi 
,=i  ’  /=i  *=i 


(2.17) 


The  displacement  vector  of  size  3x1  can  be  expressed  as  the  product  of  a 
3x24  shape  function  matrix  and  a  24x  1  nodal  displacement  vector  as  shown  in 


Appendix  C.  Taking  the  first  partial  derivatives  of  the  shape  fimctions  with  respect 
to  the  coordinates  results  in  the  following  set  of  expressions : 
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(2.18) 


dx  .  ,  dx 

i=l 

du  _ ■y  8Hi 

dy  “  dy 

du  _ ^  SHi  ^ 

dz  ,  dz 

i=\ 


dHi 

dw 

^  dHi 

’  dx 

=^,  If”'' 

dHi 

dw 

-yV/ 

dy 

=5 

dHi 

dw 

IT^' 

’  dz 

=  .2  — W' 

/=1 

The  product  of  the  6x3  partial  differential  matrix  of  Equation  (2.16)  and 
the  3x24  shape  function  matrix  are  typically  combined  into  one  6x24  matrix. 

This  matrix  is  referred  to  as  the  "B"  matrix  in  this  development.  The  "B"  matrix 


can  be  partitioned  into  sub-matrices  labeled  Bj,  where  i— 1  to  8. 

[B]  =  [[Bi][B2][B3][B,][B5][B6][B,][Bs]] 


where 


^00 

OX 

0  f  0 

0  0 

dHi  dHi  P) 
dy  dx  ^ 

n  M  M 

^  dz  dy 
dHi  rv  dHi 
dz  ^  dx 


where  z  =  1  to  8 


(2.19) 


(2.20) 


The  Galerkin  method  takes  the  shape  functions  as  the  weighting  functions 


=  W2  =  W3  = 


(2.21) 
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Actually,  the  above  equation  can  be  expressed  in  simplified  form  as 

[Ke]{d}  =  {Fe}  (2.24) 

For  a  complex  geometry,  the  finite  element  method  employs  a  mapping 
between  two  different  coordinate  systems.  In  order  to  map  a  regular  shaped 
element  in  the  'rst'  coordinate  onto  an  irregular  shaped  element  in  the  'xyz' 
coordinate,  the  Jacobian  matrix  is  required.  Equation  (2.25)  shows  the  Jacobian 
for  mapping  between  'rst'  system  and  'xyz'  coordinate  system. 

obc  ^  gz 
dr  dr  dr 

r/i=  ^  ^  fe  (2.25) 

^  ^  fe 

dt  dt  dt 
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In  terms  of  the  shape  functions  and  nodal  points,  the  Jacobian  is  expressed  as 


[J\ 


dHi 


^  uri  i 


;=1 


/=l 

8 


dr 

dH, 
ds  ' 


Y.&X- 

h  ' 


dHi 


8 

I 

/=! 

|.  dH, 

/=1 


8 


v-  dHi 

2  —yi 

i=\ 


8 

I 
(=1 
8 
y 
/=! 

h  ' 


dHi^ 

dr" 

djh. 
ds ' 


(2.26) 


Letting  [J]”'  =  [F],  the  derivatives  of  shape  functions  with  respect  to  the 
’xyz’  system  are  now  expressed  in  Equation  (2.27)  with  respect  to  the  'rst' 

coordinate  system. 


dHi  T-  fJni  -p  uiii  p 

-^  =  rii— +ri2-^+ri3  3, 


dHi 


dHi 


dHi 


dHi  p  om  p  uni  p  uii, 

-^  =  r2i-^+r22— +  r23— 


^  r  ^ 

dr  '  as  ^  23  Qt 


dHi 


(2.27) 


^  =  r3i-^+r32-^+i  33-^ 


dHi 


dHi 


dHi 


dz 


dr 


ds 


The  integral  in  terms  of  the  'xyz'  coordinate  can  be  transformed  into  the 
integral  in  terms  of  the  'rst'  coordinate  using  the  Jacobian  matrix. 

Iv  {BV{D\{BWix,y,z){d}  =  JIf,  {Bf{D\{B\\J\dV{r, s, t){d}  (2.28) 
where  |J1  is  the  determinant  of  the  Jacobian  matrix. 

The  transformation  to  the  'rst'  coordinate  system  results  in  a  sunplified 
integral  because  the  element  domain  is  regular  in  the  'rst'  coordinate  system.  The 
resultant  elements  are  called  isoparametric  elements.  The  term  isoparametric  refers 
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to  the  equality  between  the  degree  of  the  equation  for  transforming  'xyz'  into  'rsf 
coordinates  and  the  degree  of  the  shape  functions  for  estimating  displacement.  The 
shape  functions  defined  in  terms  of  the  'rsf  coordinate  for  the  eight-noded  brick 
element  are  shown  in  Equation  (2.29) 

i/i=|(l+r)(l+5)(l+0 

//2  =  i(l+r)(l-5)(l+0 

//3=i(l +/•)(! -5)(l-0 

//4  =  f(l+r)(l+5)(l-0 

//5  =  |(1-/-)(1+5)(1+0  (2.29) 

^8  =  i(l-A-)(l+5)(l-0 

The  Gauss  quadrature  is  used  for  numerical  integration.  The  volume 
integral  is  replaced  by  the  triple  summation  over  the  number  of  integration  points 
(NIP)  of  the  integral  evaluated  at  Gauss  integration  points  (r,  s,  t)  and  multiplied 
by  weighting  factors  as  shown  in  Equation  (2.30) 

NIPNIPNIP  rr 

\L{BVlDm\JW(r,sMd)  =2  E  S  [Bf[Dm\AW.WiWu{d)  (2.30) 

•'"V  i-\  i=\ 

The  results  of  the  numerical  integration  may  vary  over  the  elements  in  the 

domain  of  the  model.  Each  of  these  results  can  be  expressed  as  a  24x24  matrix 

which  is  termed  the  elemental  stiffiiess  matrix  [Ke].  The  surface  traction 
boundary  conditions  are  from  the  right  hand  side  of  Equation  (2.23).  The 
integration  of  the  applied  traction  over  the  element  surface  area  results  in  the 
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24  X 1  element  force  vector  {Fe}  which  represents  the  lumped  forces  at  the  nodal 


forces  in  the  x,  y,  z  directions.  Once  the  element  matrics  and  vectors  are 
determined,  they  are  assembled  into  the  global  system  of  equations 

[A^54^54  (c/t/)  54x1  =  (FF)  54x1  (2.31) 

In  order  to  solve  the  matrix  Equation  (2.31),  the  essential  boundary 
conditions  must  be  applied  in  advance.  In  this  study,  since  the  cell  is  symmetric, 
symmetric  boundary  conditions  are  applied  to  the  symmetric  planes.  In  addition, 
for  the  displacement  boundary  control,  as  shown  in  Fig  4,  a  uniform 
displacement  is  applied  in  the  micromechanics  model.  For  this  case,  the  elastic 
modulus  can  be  calculated  from 

E  =  -^  (2.32) 

where  Gc,  composite  average  stress,  can  be  calculated  from  subcell  stresses  and  Sc, 
composite  average  strain,  calculated  from  the  given  displacement  condition.  For 
another  case,  called  force  boundary  control,  a  uniform  traction  is  applied  instead 
of  the  displacement.  For  this  condition,  Oe  can  be  obtained  from  the  given 
force.  Then  from  the  displacement  change  of  each  node  of  the  subcell,  we  can 

calculate  the  subcell  strain.  The  cell  strain  Sc  is  obtained  by  averaging  the  subcell 
strains. 
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(b) 


Figure  4  (a)  Displacement  boundary  control  applied  on  longitudinal  or 
transverse  direction 

(b)  Force  boundary  control  applied  on  longitudinal  or  transverse 
direction 
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2.  Interface  Crack 


If  there  is  a  partial  crack  between  interface  of  the  fiber  and  matrix,  we 
assume  there  is  a  spring  between  the  connected  nodes  (between  the  fiber  and 
matrix  elements)  in  each  'u',  V,  V'  direction.  The  spring  constants  are  called 


kx,  ky,kz.  The  stif&iess  matrix  for  the  spring  element  is 

^  kx  0  0  -kx  0  o' 

0  it;,  0  0  -ky  0 

0  0  A:;,  0  0  -kz 

-kx  0  0  kx  0  0 

0  -kyO  0  ky  0 

^  0  0  -itz  0  0  kz^ 


(2.33) 


We  add  this  stif&iess  matrix  into  the  global  matrix.  This  matrix  equation 
can  calculate  the  stress  and  strain  of  the  composite  and  obtain  the  composite 
material  property  for  partial  interface  crack  between  the  fiber  and  matrix. 
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C.  MATLAB  COMPUTER  PROGRAM  DEVELOPMENT 


The  previously  derived  finite  element  formulation  was  programmed  using 
the  MATLAB  program  to  simulate  the  composite  material  property  for  perfect 
bonding  and  a  partial  interface  crack  between  the  fiber  and  matrix.  The  main 
procedure  of  the  program  is  as  follows ; 

1.  Input  X,  y,  z  coordinates  and  connectivity  to  denote  the  area  of  the  fiber  and 
matrix  using  NODAL. 

2.  Input  the  material  properties  of  fiber  and  matrix,  and  calculate  the  subcell 
material  property  matrix  [D]  using  MATER. 

3.  LFsing  the  Gauss-integration  technique  to  calculate  the  Jacobian  matrix, 
obtain  [B]’^  and  [B]  in  the  'rst'  coordinate  system  using  GAUSS  and 
JACOBIAN. 

4.  Compute  the  local  stiffiiess  matrix  [KJ  from 

5.  Assemble  the  local  stif&iess  matrix  [KJ  of  each  subcell  into  the  system 
stiffriess  matrix  [KK]  using  GLOBKF  Action. 

6.  Apply  boundary  conditions  to  modify  the  global  stiffiiess  matrix  using 
MOD. 

7.  For  the  interface  crack  analysis,  the  crack  stiffiiess  matrix  is  added  to  the 
system  stiffiiess  matrix  before  applying  the  boundary  condition. 

8.  Post-process  to  calculate  stress  and  strain  in  each  subcell  under  either 
displacement  or  force  boundary  control,  and  calculate  the  composite  stress 
and  strain  to  obtain  the  smeared  material  properties. 
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III.  RESULTS  AND  DISCUSSION 


Graphite/epoxy  is  a  continuous  fibrous  composite  material.  Graphite  fibers 
were  considered  to  be  transversely  isotropic  and  the  epoxy  matrix  was  considered 
to  be  isotropic.  The  corresponding  properties  are  listed  in  Table  1.  The  effective 
moduli  of  the  composites  predicted  by  three  dimensional  finite  element  model 
were  compared  to  the  experimental  data  [Ref.  6]. 

For  the  three  dimensional  finite  element  method,  elastic  and  shear  moduli 
in  both  longitudinal  and  transverse  directions,  respectively,  are  compared  as  well 
as  the  transverse  Poisson's  ratio.  The  experimental  and  computed  material 
properties  of  the  graphite/epoxy  composite  show  good  agreement  between  them. 
The  comparisons  are  shown  in  Figure  5  to  Figure  9.  In  general,  either  the 
displacement  boundary  control  or  the  force  boundary  control  produces  similar 
results.  However,  there  is  a  large  difference  on  longitudinal  elastic  modulus.  In  the 
transverse  direction,  the  predicted  results  are  close  each  other.  Since,  in  the 
transverse  direction,  the  acted  plane  is  matrix  which  is  much  softer  than  the 
fiber,  the  cell  deformation  depends  on  the  matrix  subcell  more  than  the  fiber  cell. 
Otherwise,  in  the  longitudinal  direction,  the  fiber  dominates  the  deformation. 

For  the  interface  crack  analysis,  the  spring  constants  kx,  ky,  Icz  are  assumed 

five  different  values,  10^‘,  5x  10‘®,  2x  10‘",  10‘“,  5x  10^  When  the  spring  constant 
is  on  the  order  of  10“ ,  it  is  assumed  there  is  no  crack  between  the  fiber  and 
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matrix  interface.  On  the  other  hand,  when  the  spring  constants  are  of  a  very  low 
order  of  magnitude,  it  means  there  is  complete  debonding  between  the  fiber  and 
matrix.  Otherwise,  there  is  a  partial  debonding.  The  results  are  shown  in  Figures 
10-14  for  the  displacement  control  model.  The  results  show  that  if  an  interface 
crack  occurs  (meaning  that  spring  constant  decreases),  the  composite  material 
property  of  the  longitudinal  elastic  modulus  Ex  will  not  be  influenced.  Both 
transverse  elastic  and  shear  moduli  decrease  as  the  interface  crack  propagates. 
Figures  15-19  show  the  solution  for  the  force  control  case.  In  this  case,  the 
interface  cracks  affect  the  longitudinal  elastic  modulus  and  increase  the  transverse 
shear  modulus,  which  does  not  make  sense  intuitively.  As  a  result,  the 
displacement  control  should  be  used  for  the  interface  crack  model. 


24 


LONGITUDINAL  ELASTIC  MODULUS 


FIBER  VOLUME  FRACTION 


Figure  5  Longitudinal  elastic  modulus  of  a  graphite/epoxy  composite 
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TRANSVERSE  ELASTIC  MODULUS 


FIBER  VOLUME  FRACTION 

Fignra  6  Transverse  elastic  modulus  of  a  graphite/epoxy  composite 
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Figure  7  Inplane  shear  modulus  of  a  graphite/epoxy  composite 
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TRANSVERSE  POISSON's  RATIO 


FIBER  VOLUME  FRACTION 

Figure  8  Transverse  Poisson's  ratio  of  a  graphite/epoxy  composite 
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Figure  9  Transverse  shear  modulus  of  a  graphite/epoxy  composite 


FIBER  VOLUME  FRACTION 


Figure  10  Longitudinal  elastic  modulus  of  a  graphite/epoxy  composite  with  interfece  crack 
(spring  constant :  1 -»  10"  2->5xl0'‘’  3 ->2x10'®  4^10'“  5 ->5x10’) 


FIBER  VOLUME  FRACTION 


Figure  11  Transverse  elastic  modulus  of  a  graphite/epoxy  composite  with  interface  crack 
(spring  constant ;  1  — >  10"  2-+5xl0'‘’  3 -♦2x10''’  4-^10'''  5 -♦5x10’) 
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Figare  13  Transverse  Poisson's  ratio  of  a  graphite/epoxy  composite  with  interface  crack 
(spring  constant :  1 -♦  10"  2 ->5x10'"  3-»2xl0'“  4->10'"  5 ->5x10’) 


otexperiment  data 
-:disp.  cntrl.  data 
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Figure  14  Transverse  shear  modulus  of  a  graphite/epoxy  composite  with  interface  crack 
(spring  constant :  1 10"  2->5xl0'“  3 -►2x10''’  4 ->10'“  5 5x10’) 


FIBER  VOLUME  FRACTION 


Figure  15  Longitudinal  elastic  modulus  of  a  graphite/epoxy  composite  with  interface  crack 
(spring constant :  1 10*’  2-^5xl0”’  3— >•2x10*®  4-^10’®  5— »5xl0®) 
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TRANSVERSE  POISSON! 


Figure  18  Transverse  Poisson's  ratio  of  a  graphite/epoxy  composite  with  interface  crack 
(spring  constant :  1 ->  10"  2->5xl0'®  3->2xl0'“  4-^10'“  5 ->^5x10’) 


IV.  CONCLUSIONS  AND  RECOMMENDATION 


In  this  study,  a  modified  analytical  micro-mechanics  model  and  3-D 
FEM  models  were  presented  in  order  to  predict  the  effective  material  properties  of 
composites  based  on  their  constituent  material  properties.  Then,  we  can  summarize 

as  follows  : 


1.  For  the  analytical  model,  the  new  model  was  slightly  better  than  the  old 
models  in  predicting  the  transverse  elastic  modulus.  Other  properties  were 
the  same  between  the  two  models. 

2  This  research  demonstrates  that  the  3-D  FEM  programmed  by  the  MATLAB 

software  can  be  efficiently  used  to  model  the  behavior  of  fibrous  composite 
material.  Comparisons  between  the  predicted  values  and  experimental 
dataproved  the  accuracy  of  the  present  model. 

3.  For  fiber  dominated  composite  materials,  using  displacement  boundary 
control  can  provide  better  results  than  using  force  boundary  control, 
especially  for  the  longitudinal  elastic  modulus. 

4  If  there  are  more  elements,  we  might  obtain  better  results  under  force 
boundary  control.  However,  using  more  elements  is  computationally 

expensive. 

5.  The  FEM  model  can  be  expanded  to  analyze  elastic  or  non-elastic  behavior. 
We  can  also  model  interfacial  cracks  between  the  fiber  and  matrix.  For  both 
cracked  and  uncracked  composites,  comparisons  between  the  present  redicted 
values  and  experimental  data  proved  the  validity  and  accuracy  of  the  present 

model. 

6.  For  a  partial  interfacial  debonding  between  the  fiber  and  matrix,  the  proper 
spring  constants  may  be  obtained  fi-om  a  more  detailed  finite  element 

analysis  of  micromechanics. 
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APPENDIX  A 


Ict  b 

■  ^22  -  ^22 

2.  <^22  ~  ^22 

3.  a^3=a33 

4.  (^22  —  ^22 

5.  aS22  +  (1  -  ^)^22  ~  ^^22  ~  (1  “  ^)^22  =  0 

6.  i7833  -a853  +  (l-a)833  -(l-a)S33  =0 

7.  8^1  -  8?,  =0 

8.  8ii  -  8fi  =  0 

9.  8n  “ 6ii  =0 

10.  ill  =a^8ii +a(l -^7)8n +a(l -a)8ii +(1 -a)^8n 

11.  822  =a^e22+a(l-a)E22  +  a(^-^)^22  +  i'^-^f^22 

12  .  833  =i7^833  +a(l-a)853+a(l-a)833  +  (l-a)^833 
*  a  =  JVf,  subscript  "a  ,b,c,d"  means  subcell. 


APPENDIX  B 


D\\  D\2  D\2  0  0  0 

D\2  D22  D23  0  0  0 

\D\  =  D\3  D23  D33  0  0  0 

0  0  0  G,2  0  0 

0  0  0  OG23O 

0  0  0  0  OG31 

J 

where 

V  =  (1  +  V23)(l  -  V23  -  2Vi2£’2/£’i) 

V21  =  V\2E2lE\  ,  V31  =  V21  ,  Vi3  =  V31 
E3=E2  ,  V32  =  V23 

Dll  =  (\  —V22)Ei/v  ,  Di2  =  V +V23)E2/v 
Du  =  Di2 

D22  =  (^1  —  Vj2^2/^1^^2/v 

D23  —  ^^23  +  Vj2^2/^1^-£^2/v 

Beacuse  G,2  is  given,  therefore 

G^23  =  ^2/2(1  +V32)  ,  G31  =  Gi2 

*  E  j  :  longitudinal  elastic  modulus  G^^  :  longitudinal  shear  modulus 
E2  :  transverse  elastic  modulus  G^^ :  transverse  shear  modulus 
Vj2  :  longitudinal  Poisson's  ratio 
V22  :  transverse  Poisson's  ratio 
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APPENDIX  C 


H\  0  0  H2OO  H3OO  Ha  0  0  H5  0  0  He  0  0  Hi  0  0  H^OO 

OHxO  OH2O  OH3O  OHaO  OH5O  OHeO  0  Hi  0  OHgO 

0  0  Hi  0  0  H2  0  0  Hi  0  0  Ha  0  0  Hs  0  0  He  0  0  Hi  0  0  i/g 
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